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ABSTRACT 

We use stripped-down versions of three semi-analytic galaxy formation models to 
study the influence of different assumptions about gas cooling and galaxy mergers. 
By running the three models on identical sets of merger trees extracted from high- 
resolution cosmological iV-body simulations, we are able to perform both statistical 
analyses and halo-by-halo comparisons. Our study demonstrates that there is a good 
statistical agreement between the three models used here, when operating on the same 
merger trees, reflecting a general agreement in the underlying framework for semi- 
analytic models. We also show, however, that various assumptions that are commonly 
adopted to treat gas cooling and galaxy mergers can lead to significantly different 
results, at least in some regimes. In particular, we find that the different models 
adopted for gas cooling lead to similar results for mass scales comparable to that of 
our own Galaxy. Significant differences, however, arise at larger mass scales. These 
are largely (but not entirely) due to different treatments of the 'rapid cooling' regime, 
and different assumptions about the hot gas distribution. At this mass regime, the 
predicted cooling rates can differ up to about one order of magnitude, with important 
implications on the relative weight that these models give to AGN feedback in order 
to counter-act excessive gas condensation in relatively massive haloes at low redshift. 
Different assumptions in the modelling of galaxy mergers can also result in significant 
differences in the timings of mergers, with important consequences for the formation 
and evolution of massive galaxies. 

Key words: galaxies: formation ~ galaxies: evolution - galaxies: cooling flows - 
galaxies: interactions. 



1 INTRODUCTION 

Understanding how galaxies form and the physics that drive 
their evolution has been a long-standing problem in modern 
astrophysics. A number of observational tests have recently 
succeeded in determining the fundamental cosmological pa- 
rameters with uncertainties of only a few per cent, thus effec- 
tively removing a large part of the parameter space in galaxy 
formation studies. We are left, however, with the problem of 
dealing with our 'ignorance' of complex physical processes, 
that are inter-twined in an entangled network of actions, 
back-reactions, and self-regulations. 

Over the past decades, different methods have been de- 
veloped to study galaxy formation in a cosmological con- 



text. Among these, semi-analytic models have turned into 
a flexible and widely used tool to provide detailed pre- 
dictions of galaxy properties at different cosmic epochs. 
Th ese techniques find t heir seeds in the pioneering work 
by IWhite fc Reed { 19781 ). were laid out in a more detailed 
form in the early 1990s jWhite fc Frenkl Il99ll : ICold Il99ll : 
iKauffmann. White fc Guiderdonilligoi T and have been sub- 
stantially extended and refined in the last years by a number 
of different groups (for a review, see lBaughll200d ). In these 
models, the evolution of the baryonic components is mod- 
elled adopting simple yet physical and/or observationally 
motivated recipes, coupled in a set of differential equations 
that describe the variation in mass as a function of time of 
different galactic components (e.g. stars, gas, metals). 
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While it is relatively easy to compare results from differ- 
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ent model^, it is more complicated to understand the origin 
of any difference or similarity between them. This difficulty 
stems primarily from the fact that different groups adopt 
different sets of prescriptions (that are equally reasonable, 
given our poor understanding of the physics at play) and 
that, as mentioned above, the final results are given by a 
combination of these. There are, however, also a number 
of more subtle differences that are more 'technical' in na- 
ture (e.g. the particular mass definition adopted, the cool- 
ing functions used, the use of analytic or numerical merger 
trees, etc.). The precise influence on models' results of at 
least some of these details is unclear. For example, it has 
been sh own that the extended Press- Schechter (EPS) for- 
malism l|Bond et al.lll99ll : lBowerlll99ll ) does not provide an 
adequate description of the me rger trees extrac t ed directly 
from numerical simulations (B enson et ahllioOSl : ICole et al.l 
I2OO8I ). Although some of these most recent studies have pro- 
vided 'corrections' to analytic merger trees, many applica- 
tions are still carried out using the classical EPS formalism, 
and little work has been done to understand to which mea- 
sure this affects predictions of galaxy formation models. 

In this paper, we compare results from three indepen- 
dently developed semi-analytic models. Our goal is not to 
predict or reproduce any specific observation. Rather our 
aim is to understand the level of agreement between differ- 
ent semi-analytic models, with a minimum of assumptions 
or free parameters. This requires (1) that the models are im- 
plemented on identical merger trees, such that results can be 
compared on a case-by-case basis and any differences can be 
attributed to specific parametrizations of physical processes; 
and (2) that a minimum of physical processes are included 
in the models in order to avoid possible degeneracies, and to 
hopefully illuminate the effects of specific parametrizations 
or parameter choices. 

In our study, the first requirement has been satisfied by 
creating a standard set of halo merger trees extracted from 
Ai'-body simulations, and running each model on these trees. 
The second requirement is somewhat more demanding as 
modern semi-analytic models contain treatment of numer- 
ous, coupled physical processes. We have chosen to simplify 
the models as much as possible by removing all physics other 
than gas cooling and galaxy mergers. This allows us to focus 
on the influence of different assumptions typically made to 
model these two physical processes, that represent two ba- 
sic ingredients of any galaxy formation model. Using large 
samples of identical haloes, we are able to compare results 
both in a statistical fashi on and on a halo- by-halo basis. 

Previous studies fe.g. Benson et al."20 0ll :l Yoshida et al.l 
l2002l : lHellv"et a l. 200a ; .Cattanco o t al. 200^ have compared 
numerical predictions from stripped-down versions of semi- 
analytic models with those from hydrodynamical simula- 
tions, to verify whether these methods provide consistent 
predictions in the idealized case in which only gas cool- 
ing is included. The general consensus is that the cooling 
model usually employed in semi-analytic models is in good 



^ A number of galaxy catalogues have been made publicly avail- 
able by various groups; results from different versions of two of 
the models used in this study are available through a relational 
database accessible at: 

[http: / /www.mpa-garching.mpg.de/millennium/l 



agreement with hydrodynamical simulations that adopt the 
same physics. More recent studies focused on object-by- 
object comparisons, however, have highlighted a number of 
important differences that were 'hidden' in the relatively 
good agreement obtain ed by previous studies focusing on 
statistical compar isons (|Saro et al.ll201ol ). A recent study by 
IViola et al.l l|2008h . in particular, has shown that the cooling 
model implemented in MORGANA (one of the models used 
in this study) predicts cooling rates that are significantly 
larger than those predicted from their im plementation of 
the 'c lassical' cooling model, proposed by IWhite &i FrenkI 
( I991I ). In addition, Viola et al. have shown that MORGANA 
provides results that are in good agreement with those of 
controlled numerical experiments of isolated haloes, with hot 
gas in hydrostatic equilibrium. While both SPH and semi- 
analytic techniques have their own weaknesses, making it 
unclear which of the two (if either) is providing the 'cor- 
rect' answer, these results appear confusing. It is therefore 
interesting to study how the different possible assumptions 
that can be made to model the evolution of cooling gas, 
propagate on predictions from galaxy formation models. 

Modelling of galaxy mergers has not been considered 
a major concern, but different assumptions about merging 
time-scales can be made, and these may have important con- 
sequences for the inferred stellar assembly history of massive 
galaxies, including brightest cluster galaxies. In addition, re- 
cent work has shown that the classical dynamical friction 
formula usually adopted in semi-analytic models tends to 
under-estimate merging times computed from controlled nu- 
merical experiments and high-resolution c osmologica l simu- 
lation s l|Bovlan-Kolchin. Ma fc QuataertI [2OO8, Jiang et al.l 
l2008h . Results from these studies, however, have not yet 
converged on the appropriate correction(s) that should be 
applied to the classical formula. 

In this paper, we will not address the issue of what is 
the best way to model galaxy mergers or gas cooling. Instead, 
we will concentrate on the differences due to alternative im- 
plementations of these physical processes, with the aim of 
understanding their effects in a full semi-analytic model. We 
will also explore what these effects might imply for the im- 
portance of other physical processes. 

The numerical simulations and merger trees used in our 
study are described in Section [S] In Section |31 we describe 
in detail how gas cooling and galaxy mergers are treated 
in each of the three models used in this work. Section |4] 
and Section [S] present our results for two halo samples. In 
Section [6l we discuss the influence of numerical resolution 
and of different schemes for the construction of merger trees. 
Section [7] compares the different implementations of merger 
times adopted in the three models considered. Finally, we 
summarize and discuss our results, and give our conclusions 
in Section [S] 



2 THE SIMULATIONS AND THE MERGER 
TREES 

This work takes advantage of two large high-resolution 
cosmolog ical s i mulati ons: the Millennium Simulation (MS, 
l2005l). 



Springcl et al.l 



and the 



^ Millennium-II (MS-II, 

iBovlan-Kolchin et ahl l2009h . The MS follows N = 2160^ 
particles of mass 8.6 x 10* h~^MQ, within a comoving box 
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Name 


-^box 


e [h-i kpc] 


mp [/i-^Mq] 


MS 


500 


5.0 


8.61 X 10* 


MS-II 


100 


1.0 


6.89 X 


mini-MS-II 


100 


5.0 


8.61 X 10* 



Table 1. Some basic properties of the three simulations used 
in this study: the side length of the simulation box Lboxi the 
Plummer-equivalent force softening e, and the particle mass mp. 



of size 500 ^"^Mpc on a side. The MS-II follows the evo- 
lution of the same number of particles in a volume that is 
125 times smaller than for the MS (100/i~^Mpc on a side). 
The particle mass is correspondingly 125 times smaller than 
for the MS (6.9 x 10® H'^Mq), allowing haloes similar to 
those hosting our Milky Way to be resolved with hundreds of 
thousands of particles. For both simulations, the cosmolog- 
ical model is a ACDM model with = 0.25, fib = 0.045, 
h = 0.73, f^A = 0.75, n = 1, and as = 0.9. The Hubble 
constant is parameterised as Ho = 100 /ikms^^Mpc^^. In 
order to test how numerical resolution affects our results, we 
will also use the mini-Mille nnium-II simulation (mini-MSII; 
iBovlan-Kolchin et aLl bOOgj). which was run using the same 
initial conditions and volume as for the MS-II, but at the 
mass and force resolution as for the original MS (the num- 
ber of particles is therefore 432'^). The basic properties of 
the three simulations used in this study are summarised in 
Table [U 

For each simulation snapshot (64 for the MS, 68 for the 
MS-II), group catalogues were constructed using a standard 
friends-of-friends (FOF) algorithm, with a linking length 
of 0.2 in units of the mean particle separation. Each group 
was then decomposed into a s et of disjoint substr uctures 
using the algorithm SUBFIND l|Springel et aLlbOOll ). which 
iteratively determines the self-bound subunits within a FOF 
group. The most massive of these substructures is often re- 
ferred to as the main halo, while this and all other sub- 
structures are all referred to as subhaloes or substructures. 
Only subhaloes that retain at least 20 bound particles af- 
ter a gravitational unbinding procedure are considered 'gen- 
uine' subhaloes, and are use d to construct merge r history 
trees as explained in detai l in lSpringel et al. ||2005|) (see als o 
iDe Lucia fc BlaizotI I2OO7I and iBovlan-Kolchin et aU l2009f ). 
The subhalo detection limit is therefore set to 2.36 x 10^°, 
and to 1.89 x 10* M© for the MS and the MS-II, respectively. 
Note that some FOF haloes do not contain 20 self-bound 
particles; such FOF haloes are not included in the merger 
trees. 

The comparison discussed below will focus on two sam- 
ples of haloes. The first sample consists of 100 haloes, se- 
lected from the MS-II on the basis of their mass at 2 = 0, 
with logM2oo between 11.5 and 12.5. Here M200 is in units 
of /i~^M0, and is defined as the mass within a sphere of 
density 200 times the critical density of the Universe at 
the corresponding redshif10. We will refer to this sample 



^ The corresponding radius, -R20O1 has been shown to approxi- 
mately demarcate the inner regions of haloes which are in dynam- 
ical equilibrium, fro m the outer regions which are still infalling 
llCole fc Lacevlll996l'l . In the following, we will refer to this ra- 
dius as the 'virial radius', -Rvir- The virial mass Myir is the mass 



as the 'Milky- Way like' sample. A second sample of other 
100 haloes was selected from the MS by taking haloes that 
have a number density of 10~^ at z ~ 2, and that end up in 
massive groups/clusters at 2: = 0. The adopted number den- 
sity is comparable to t hat of submillimiter galaxies at z ~ 2 
l|Chapman et al.ll2004l . and references therein). We will refer 
to this sample as the 'SCUBA-like' sample. 

As mentioned above, trees for the MS and MS-II were 
originally constructed at the subhalo level. Many semi- 
analytic models, however, are based on FOF merger trees. 
We have therefore constructed FOF-based merger trees for 
our halo samples. There are many possible ways to construct 
FOF merger trees from subhalo trees. We have chosen one of 
the most straightforward methods: since there exists a one- 
to-one correspondence between FOF haloes in the subhalo 
trees and main subhaloes, we are free to assign the merger 
history of a FOF halo to the merger history of its dominant 
subhalo. We use this mapping from dominant subhaloes to 
FOF haloes to directly link the merger history of dominant 
subhaloes to the merger history of their host FOF haloes. 
So, for example, if dominant subhalo a of FOF group A 
has the descendant j3 in FOF group B, then FOF halo B 
is defined to be the descendant of FOF halo A in our FOF 
trees. 

More sophisticated algorithms for constructing FOF 
trees from the MS or MS-II s ubhalo trees exist (e.g. 
iFakhouri fc Mall2008l : [Genel et al.|[2 008'). In particular, these 
algorithms filter out "unphysical" FOF mergers - for exam- 
ple, those due to chance associations of two FOF haloes for 
only one snapshot. Each of these algorithms has different 
strengths and weaknesses. The primary virtue of the algo- 
rithm we have chosen is its simplicity, and the advantage 
that each FOF tree can be easily connected to the corre- 
sponding subhalo tree (we will use subhalo-based merger 
trees in Section [6]). 



3 THE MODELS 

In this study, we use three different and independently de- 
veloped semi-analytic codes. In the following, we will refer to 
them as the Munich model, the Durham model, and MOR- 
GANA. As explained above, we have used stripped-down ver- 
sions of the models that only include gas cooling and galaxy 
mergers, so as to focus on a few specific aspects of the mod- 
elling. In addition, all models have been adapted to run on 
the merger trees described in Section (2] In the following, we 
describe in more details how gas cooling and galaxy mergers 
are treated in each of the models used in our study, and the 
changes that were made in order to adapt each code to the 
same merger trees. 



3.1 The Munich model 

The version of t he Munich model used in this study is the 
one described in lDe Lucia fc BlaizotI (|2003), and we refer to 
the original paper and references therein for more details. 
The rate of gas cooling is computed following the model 



contained within the sphere defined by this radius, and the virial 
velocity Wir is the circular velocity at iJvir. 
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originally proposed IWhite fc Freri^ |l99lf). and an imple - 
mentation similar to that adopted in Springel et alj l|200lh . 
The Munich model assumes that the hot gas within dark 
matter haloes follows an isothermal profile: 

, _ Mhot 

For each new snapshot, the total amount of hot gas 
available for cooling in each halo is estimated as follows: 



/bMv 



/ J CO 



) 

old 



(1) 



where the sum extends over all the galaxies in the FOF halo 
and /b is the baryon fraction of the Universe, for which we 
assume the value 0.01/0. Eq.[T] can provide, in a few cases, a 
negative number (this occurs typically after important halo 
mergers). In this case, the amount of hot gas is set to zero, 
and no cooling is allowed in the remnant halo. 

The equations driving galaxy evolution are then solved 
using 20 time-steps between each pair of simulation snap- 
shots. A local cooling time is defined as the ratio between 
the specific thermal energy content of the gas and the cool- 
ing rate per unit volume: 



2fim-pnl{r)K{T, Z) 



(2) 



In the above equation, fim^ is the mean particle mass, ndr) 
is the electron density, k is the Boltzmann constant, and 
A(r, Z) represents the cooling rate. The latter is strongly 
dependent on the virial temperature of the halo, and on the 
metallicity of the gas. In the Munich model, these depen- 
dencies are accou nted for by using the c oUisional ionization 
cooling curves by [Sutherland fc Dopital (|l993i' ). Since chem- 
ical enrichment is switched off in this study, we only use the 
calculation corresponding to 'primordial' composition. The 
virial temperature of the halo is determined using the hydro- 
static equilibrium equation, and relates the gas temperature 
to the circular velocity of the halo: 



1 AtmH..2 

2 k 



or Tvir = 35.9(Vvir/kms~^)^ K 



where mn is the mass of the hydrogen atom, and /i is the 
mean molecular mass. 

A 'cooling radius' is then computed as the radius at 
which the local cooling time is equal to the halo dynamical 
time. We note that in the original work by White & Frenk, 
the cooling radius was defined equating the local cooling 
time to the age of the Universe, which is about one order 
of magni tude larger than the halo dynamical t ime. As dis- 
cussed in lDe Lucia. Kauffmann fc Whit^ (|2004l ). the partic- 
ular choice currently adopted in the Munich model was mo- 
tivated by the significant enhancement of cooling rates when 
adopting metal dependent cooling functions. 

If the cooling radius lies within the virial radius of the 
halo under consideration, the gas is assumed to cool quasi- 
statically, and the cooling rate is modelled by a simple inflow 
equation: 

dMcooi 2 drcooi 

= 47rpg(rcooi)r-cooi- 



dt 



At 



^ In this work, all models assume the same value for the universal 
baryon fraction. 



At early times, and for low-mass haloes, the formal cool- 
ing radius can be much larger than the virial radius. In this 
case, the infalling gas is never expected to come to hydro- 
static equilibrium, and the supply of cold gas for star for- 
mation is limited by the infall rate. In this 'rapid cooling 
regime', we assume that all new diffuse gas that is accreted 
onto the halo is immediately made available for star forma- 
tion in the central galaxy of the halo under consideration. 

In its standard implementation, the Munich model fol- 
lows explicitly dark matter haloes when they are accreted 
onto larger systems. This allows the dynamics of satellite 
galaxies residing in infalling structures to be properly fol- 
lowed, until their parent dark matter substructures are com- 
pletely destroyed by tidal truncation and stripping (e.g. 
iDe Lucia et al.ll2004l ). When this happens, galaxies are as- 
signed a residual surviving time that is estimated from the 
relative orbit of the two merging objects, at the time of sub- 
halo disruption, using the following implementation of the 
Chandrasekhar dynamical friction formula: 



1.17 D-" M„ 



fudge 



InAdt K 



Ms: 



■ '''dyn 



(3) 



In the above equation, D is the distance between the merg- 
ing halo and the centre of the structure on which it is ac- 
creted, _Rvir is the virial radius of the accreting halo, Msax is 
the (dark matter) mass associated with the merging satel- 
lite, and Mniain is the (dark matter) mass of the accreting 
halo. The dynamical time of the halo, rdyn, is given by 

1/2 



''"dyn 



-Rvir 
Vvir 



n3 
^vir 

GMvir 



(4) 



Note that with the definition of virial mass adopted here, 
Tdyn ~ 0.1/ H(z) aud is independent of the halo mass. The 
Coulomb logarithm Adf is taken to be 1 -I- Mmain/Msat. 

For the purposes of this analysis, we have used Eq. [3] 
also for FOF-based merger trees, without adding any addi- 
ti onal orbital depende n ce, an d using _Rvir in place of D. As 
in lDe Lucia fc BlaizotI (|2007f ). we have assumed /fudge ~ 2. 
This was originally introduced to reduce the slight excess 
of bright galaxies otherwise pr oduced, and w a s mot ivated 
by some preliminary work by ISpringel et al.l (|200ll ). who 
noted that merging times inferred from Eq. [3] are typi- 
cally shorter than those directly measured using higher res- 
olution numerical simulations. As mentioned in Section [T] 
more recent work has confirmed that the classical dynami- 
cal friction formula tends to under-estimate m erging times 
iBovlan-Kolchin et al.|[2008l : I Jiang et al.ll2008l '). While these 
studies have not yet converged on the proper adjustment(s) 
to the classical dynamical friction formulation, they suggest 
that the appropriate correction(s) cannot be absorbed in 
a fudge factor in fron t of Eq. El For example, while the 
fudge factor adopted in lDe Lucia fc Blaizot (|2007|) is in good 
agreement with findings from Bovlan-Kolchin et al.l (|2008l ) 
for A/sat /A/host ~ 0.1, it is significantly lower than what 
found for smaller mass ratios. 



3.2 The Durham model 

The versi on of the Durham m odel used in this study is de- 
scribed in lBower et al.l (|2006f ). For full details on the cool- 
ing and merging times modelling , th e reade r is referred to 
iBenson et al.. (,200a ') and Cole etlo] (|2000l l. A few small 
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modifications were necessary in order to make the model 
run on the merger trees described in Section [2] These are 
briefly discussed below. 

The hot gaseous component in dark matter haloes is 
assumed to have a density profile described by the /3-model: 



Po 



[1 + (r/reo.e)2)]3/3/2 



where po is the density at the centre of the halo, Tcore is the 
radius of the core, and /3 is a parameter that sets the slope 
of the profile at radii larger than rcorc- The model assumes 
a gas density profile with rcore ~ 0.07 • -Rvir and = 2/3 
for a ll haloes, in absence of energy injection (|Benson et al.l 
I2OO3I) . Since feedback is switched off in this study, these 
are the model parameters adopted for our comparison. The 
temperature profile of the gas is assumed to be isothermal 
at the virial temperature. 

A set of 'halo formation events' are defined throughout 
each merger tree, and cooling calculations are begun and 
reset at these formation events. In particular, each halo with 
no progenitor is fiagged as a 'formation event'. For all other 
haloes in the tree, a 'formation event' fiag is assigned to all 
those for which their mass is equal to, or larger than, twice 
the mass of the halo at the previous formation event in that 
branch. For these haloes then, formation events correspond 
to halo mass doublings. 

At each new snapshot, the mass of gas that falls onto 
the halo is given by: 



Minfaii = max 



(Mvi.-^M;i,)./b,0.0 



where the sum extends over all the progenitors of the halo 
under consideration. Although Minfaii is accumulated for 
each snapshot, it is only added to the hot gas of the halo at 
a formation event. 

As for the Munich model, the following calculations are 
then done using 20 time-steps between each pair of simula- 
tion snapshots. A local cooling time is defined using Eq. [21 
and t he metal free coolin g function computed from CLOUDY 
v8.0 jFerland et al.lll998l ). A cooling radius is then computed 
by equating the local cooling time to the time since the last 
formation event (the cooling radius is not allowed to exceed 
the virial radius). An infall radius is also computed, as the 
minimum between the cooling radius and the free-fall radius 
(i.e. the radius within which gas has had time to fall ballis- 
tically to the halo centre, assuming that it began at rest at 
the previous formation event). The mass of gas that infalls 
onto the central galaxy during each time-step is finally com- 
puted as the difference between the hot gas enclosed within 
the current infall radius and the mass that was inside the 
infall radius at the previous time-step. 

Galaxy mergers are treated in a way similar to that 
done in the Munich model, with a few differences. When a 
new halo forms, each satellite is assumed to enter the halo 
on a random orbit (the most massive pre-existing galaxy 
becomes the central galaxy of the remnant halo). The dy- 
namical fr iction formulation adopted in the Durham model 
IS given m ICole et al.l (|2000l ) and reads as: 



fudge G^orbit 



0.3722 M„ 



In Adt 



■ '''dyn 



(5) 



orbits, and A/sat is the mass of the satellite galaxy includ- 
ing the mass of the dark matter halo in which it formed. 
The Coulomb logarithm Adf is taken to be Mhaio/A^aat. The 
orbital dependence is contained in Gorbit, modelled as a 
log-normal distribution with mean < log Gorbit >= —0.14 
and dispersi o n < ( logSprbit— < logSorbit >)^ >i/2— o.26. 
[Bower et ajj (|2006l ). and the stripped-down version used in 
this study, assume the value 1.5 for the dimensionless param- 
eter /fudge- Merger times are reset at each formation event, 
re-extracting orbital parameters for each satellite. 

3.3 MORGANA 

All details about the modelling of gas cooling and 
galajcy mergers adopted in MORG ANA can be found in 
|m onaco. Fontanot fc Taffonil (|2007l ). and we refer to this pa- 
per for full details. We note that in its original formulation, 
MORGANA us es merger trees obtained using PINOCCHIO 
([Monaco et al.l 2002). This algorithm, based on Lagrangian 
perturbation theory, has been shown to provide mass as- 
sembly histories of dark matter haloes that are in excellent 
agree ment with results from numerical simulations ([Li et al.[ 
[2OO7D . For the purposes of this study, we have adapted MOR- 
GANA to run on the numerical merger trees described in 
Section [2I This required some small modifications that are 
described below. We note that PINOCCHIO does not provide 
information on dark matter substructures, so MORGANA is 
essentially based on FOF merger trees. 

The hot halo phase is assumed to be spherical, in hydro- 
static equilibrium in a NFW halo, described by a polytropic 
equation of state with index 7^ = 1.15, and is assumed to fill 
the volume between the cooling radius and the virial radius 
of the halo, where it is pressure balanced. The equilibrium 
configuration of the hot halo gas is computed at each time- 
step, assuming that the gas re-adjusts quasi-statically to the 
new equilibrium configuration, in absence of major mergers. 
Under these assumptions, one obtains: 



Pgir) = PgO 1 - a 1 



ln(l + c„f„.T)\V'''^^"'' 



rg(r) = r,o i-a 1- 



ln(l + Cnfwa:) 

Cnfw 



where 



3(7p 



<^nfw(l ~l~ <^nfw ) 



VOlp V (1 + Cnfw) ln(l + Cntw) - Cnfw 

In the above equations, Cnfw ~ rhaio/rs, and x = r/rs, where 
rs is the scale radius of the NFW halo. The constants pgo and 
Tgo are defined as the extrapolations to r = of the density 
and temperature profile, while rjo is the extrapolation to 
r = of the function rj{r) — Tg{r)/Tvir- The halo virial 
temperature T^ir is defined as l/SpmnV^i^/k. 

At each new snapshot, the mass of gas that falls onto 
the halo is computed as: 



Minfaii = fh ■ max 



AU 



(^Ar + Mv7r),o.o 



where A^haio is the mass of the halo in which the satellite 



where the sum extends over all progenitors of the halo that 
are not its main progenitor, and M^^^ is the maximum virial 
mass of the main progenitor, considering all previous snap- 
shots. The infall rate of new gas is assumed to be constant 
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over the time interval between each pair of simulation snap- 
shots. The equations below are then solved using a Runge- 
Kutta integrator with adaptive time-steps. 

The cooling rate of a shell of gas of width Ar, at a 
radius r, is computed as: 



rfMcooi(r) 
dt 



Aixr 



(r)Ar 



tcool(r) 

where fcooi is the local cooling time, computed as in Eq.[2] As 
in the Munich model, cooling rates ar e computed using the 
collisi onal ionization cooling curves bv lSutherland fc Dopital 
l| 19931 ). with primordial composition. The mass deposition 
rate is then computed by summing up the contributions 
from all mass shells. The summation is carried our by taking 
into account the radial dependence of the gas density, and 
provides the following total mass deposition rate: 



dt 



^cool.O 



1 - a 1 - 



ln{l+t) 



t 



n 2/(7p-i) 



where tcooi,0 is computed using the central density pgo and 
the temperature of the gas at rcooi- 

Finally, by equating the mass cooled in a time interval 
dt with the mass contained in a shell dr, one obtains the 
evolution of the cooling radius, that is treated as a dynamical 
variable in this model: 



dr^. 



dt 



dM^ooi/dt 
47rpg(rcooi)r2„oi 



where Cs is the sound speed computed at rcooi, and is added 
to the right hand side of the above equation to allow the 
'cooling hole' to close at the sound speed. 

The cooling calculation is started when a halo appears 
for the first time, with rcooi = 0.001 • r^, and is reset after 
any halo major merger (Msat/Mmain > 0.2). Finally, cooled 
gas is incorporated onto the central galaxy at the following 
rate: 



dt 



(rcooi 



where ridyn is treated as a free parameter (we assume a value 
of 0.3 for this parameter, as in the standard MORGANA 
model) and represents the number of dynamical times, com- 
puted at the cooling radius, required for the cooled gas to 
be incorporated onto the central galaxy. In the comparison 
discussed below, the above delay in the incorporation of the 
cooled gas is neglected. 

When a halo is accreted on a larger structure, orbital 
parameters are extracted randomly from suitable distribu- 
tions that are based on results from numerical simulations. 
In particular, the eccentricity of the orbit (e — J/Jc, where 
J is the initial angular momentum of the orbit and Jc is 
the angular momentum of a circular orbit with the same en- 
ergy) is extracted from a Gaussian with mean 0.7 and vari- 
ance 0.2, while the energy of the orbit {xc = rc/rh, where 
rji is the halo radius and Vc is the radius of a circular orbit 
with the same energy) is assumed to be 0.5 for all orbits. To 
model galaxy mergers, [Monaco et al. I (l2007l) use aslight up- 
date of the formulae provided bv lTaffoni et all (|2003l ). that 
take into account dynamical friction, mass loss by tidal strip- 
ping, tidal disruption of subhaloes, and tidal shocks. Galaxy 



merger times are computed interpolating between the case 
for a 'live satellite' (the object is subject to significant mass 
losses) and that of a 'rigid' satellite (no mass loss): 



'merge, live 



(Cl,/: 



0.12 I f /•2 \ 



"^dyn 
/"sat 

X (0.25 f-it + 0.07/nfw + 1.123) (0.4 + ^-^e - 0.2)) 

r,ncrge,rigid = 0.46 • (1.7265 + 0.0416 Cnfw) (6) 

/sat In Adf 

In the above equations, /sat = A/sat/Mvir and /nfw ~ 
Csat/Cnfw Afgat and Csat are the corresponding virial quan- 
tities for the satellite halo, and are polynomial functions 
of Xc, whose e xpres sions can be found in appendix A of 
[Monaco et al.l l|2007h . The Coulomb logarithm is given by 
Adf = 1 -|- l//sat. Merger times are reset after each halo ma- 
jor merger, re-extracting orbital parameters for each satellite 
galaxy. 

3.4 Model differences and similarities 

The previous sections illustrate that the implementations of 
gas cooling and galaxy mergers in the three models used in 
this study differ in a number of details. Both the Durham 
and the Munich models adopt variation s of th e original cool- 
ing model proposed bv lWhite fc Fren^ l|l99ll ). but they use 
different gas profiles, and different definitions of the 'cooling 
time' (used to calculate the cooling radius). In addition, the 
cooling calculation is reset in the Durham model at each 
'formation event', and the Munich model assumes very ef- 
ficient cooling in the rapid cooling regime. Finally, cosmo- 
logical infall of gas onto the halo occurs at a constant rate 
between each pair of snapshots in the Munich model, while 
occurs at formation events for the Durham model. A cored 
density profile, like the one adopted in the Durham model, 
is expected to give lower cooling rates than an isothermal 
profile. It is not clear, however, how this naive expectation is 
affected by the other different assumptions discussed above. 

The adopted modelling for merger times is also very 
similar in the Durham and Munich models, as well as for 
the 'rigid' case in the MORGANA model. Indeed, by com- 
paring Eqs.OO and [S] with D — i^vir in the Munich model, 
©orbit = 1 in the Durham model, and cnfw ~ 10 in (the 
rigid version of) MORGANA, the implementations differ only 
in the numerical pre-factor (Munich: 1.17 /fudge; Durham: 
0.3722 /fudge; MORGANA: 0.348) and in the implementation 
of the Coulomb logarithm. (The full 'live' version of dynam- 
ical friction in MORGANA is somewhat more complicated.) 

MORGANA differs significantly from the Durham and 
the Munich models, both in its gas cooling and merger times 
implementations. As explained above, in the Durham and 
Munich model the cooling radius is computed by equating 
the local cooling time (Eq. [2]) to some given time (the halo 
dynamical time in the Munich model and the time since 
the last formation event in the Durham model). MORGANA 
computes the cooling rate at each mass shell, integrates over 
the contribution of all shells, and follows the evolution of the 
cooling radius by assuming that the transition from the hot 
to the cold phase is fast enough to create a sharp edge in 
the density profile of the hot gas. The cooling radius 'closes' 
at the sound speed. 

As mentioned in Section [TJ IViola et all Hooi) have 
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shown that results from this model are in good agree- 
ment with controlled hydrodynamical simulations of isolated 
haloes, with hot gas in hydrostatic equilibrium in a NFW 
halo. Viola et al. have also shown that their implementation 
of the classical model underpredicts the amount of gas cool- 
ing with respect to the model adopted in MORGANA, par- 
ticularly at early times. Their implementation of the 'clas- 
sical' model differs in detail from those adopted in the Mu- 
nich and Durham models, however. In addition, it is un- 
clear that results obtained from their controlled simulations 
should remain valid when using cosmologically motivated 
halo merger histories, as done in this study. The formulation 
adopted to model dynamical friction in MORGANA is based 
on controlled simulations and analytic models, and takes 
into account dynamical friction, mass loss by tidal strip- 
ping, tidal disruption of subhaloes, and tidal shocks. Since 
the other two models treat satellites as rigid systems, they 
are expected to provide somewhat shorter merger times. We 
will come back to a more detailed comparison between the 
adopted formulations in Section [T] 

MORGANA and t he Munich model compute cooling 
rates using results by [Sutherland fc Dopital (|l993l ). while 
the Durham model adopts updated rates computed us- 
ing CLOUDY. We have verified that the metal free cool- 
ing function computed using CLOUDY does not differ 
significantly from the prim ordial cooling function from 
[Sutherland fc Dopital l|l993r) . It should be noted, however, 
that this is not the case for non-zero metallicities, where the 
Sutherland & Dopita calculation tends to over-estimate the 
cooling rates computed using CLOUDY. 

Finally, we note that in all three models used in this 
study, gas cooling occurs only on central galaxies. When 
galaxies are accreted onto a more massive system, their hot 
reservoir is assumed to be instantaneously stripped and as- 
sociated with the parent halo of the remnant central galaxy. 



4 MILKY- WAY HALOES 

In this section we will analyse results of the three mod- 
els described above for a sample of Milky- Way like haloes. 
As discussed in Section [21 these haloes have been selected 
only on the basis of their present virial mass, with no addi- 
tional constraint on their merger history or isolation. An in- 
depth analysis of the full sample of Milky Way-mass haloes 
in the MS-II (comprised of o ver 7600 haloes) is presented in 
iBovlan-Kolchin et al.l (|2010l '). 

Fig. [T] shows the dark matter mass scaled by the uni- 
versal baryon fraction (top panel), the evolution of the cold 
gas (second panels from top) and hot gas (third panels from 
top) components associated with the central galaxy, and the 
cooling rate (bottom panels) for two MW-like haloes. The 
cold and hot components have been normalized by the dark 
matter mass of the parent halo, which is the same in all three 
models by construction. Red, green and blue lines show re- 
sults from the Durham, Munich and MORGANA models, re- 
spectively. The haloes chosen for this figure provide two rep- 
resentative examples for a MW halo with a rather quiet mass 
accretion history (left panels) and one with a larger number 
of massive mergers (right panels). To avoid complications 
arising from a different treatment of merging times, the evo- 
lution of the cold gas content is shown for both tmrg = oo 



(solid lines) and tmrg ~ (dashed lines). The amount of hot 
gas associated with the central galaxy, as well as the cool- 
ing rate, are not affected by the particular merging model 
adopted, because the hot gas associated with galaxies falling 
onto a larger system (i.e. becoming satellites) is instanta- 
neously stripped and associated with the hot gas component 
of the main halo, in all three models. 

Fig. [T] shows a quite good degree of agreement between 
the three models used in this study. At high redshift, when 
the halo is in a rapid accretion regime, a large fraction of 
its baryonic mass is converted into cold gas in the Munich 
model and in MORGANA, while cooling appears to be less 
efficient in this regime in the Durham model. By redshift 
~ 7 — 9, all models converge to about the same cold gas and 
hot gas fractions. At lower redshift, the evolution of both 
components is very similar, with small differences between 
different models at present (see below). As expected, the 
evolution of both baryonic components is more noisy for 
the halo whose mass accretion history is characterized by 
a larger number of important mergers. In particular, the 
Durham model shows a quite noisy behaviour in the hot 
gas evolution, which is due to the fact that gas infalling 
onto the halo is added to the hot gas component only at 
'formation events' (see Section [S}. For the halo shown in the 
right panels, the differences between the amounts of cold and 
hot gas predicted by the three models are somewhat larger 
than for the halo with a smoother accretion history shown 
in the left panels, and differences appear to accumulate at 
each merger event. 

The predicted cooling rates are very noisy for all three 
models used in this study, and in both examples shown 
in Fig. [T] In these, the highest cooling rates are obtained 
in the Durham model and are of the order of ~ 1.8 x 
lO^'^ M0Gyr~^. In a number of the other haloes in the MW 
sample, cooling rates as high as about twice this value are 
obtained, and there is some tendency for the Durham model 
to provide the noisiest behaviour, in particular at very low 
redshift {z < 0.5). Overall, however, the evolution of the 
cooling rates predicted by the three models is quite similar. 

In order to study any systematics in the evolution pre- 
dicted by the three models used in this study, we have run 
each code on the total sample of 100 MW-like haloes. Fig. [5] 
shows the average evolution of the cold (left panel) and hot 
(middle panel) gas components associated with the central 
galaxy, and of the cooling rate (right panel). In the left 
panel, solid and dashed lines correspond to the tmrg = oo 
and tmrg = run, respectively. Red, green and blue lines 
show again results from the Durham, Munich and MOR- 
GANA models, respectively. The fraction of cold gas associ- 
ated with the central galaxy rises steeply in all models and 
at progressively decreasing redshift for the Munich model, 
MORGANA, and for the Durham model. The later evolution 
of this component is very similar in all three models when 
merger times are set to zero, with a slight tendency for the 
Munich model to predict the largest cold gas fraction and 
the lowest hot gas fraction at present (see below). When 
merger times are set to oo, the Munich model predicts the 
highest cold gas fractions at early times and the lowest cold 
gas fractions at low redshift. MORGANA has a similar be- 
haviour but the cold gas fraction rises slightly later than in 
the Munich model at early times, and is slightly larger than 
predicted by the Munich model at late times. Finally, in the 
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Figure 1. From top to bottom; dark matter, cold gas, hot gas, and gas cooling rate for two MW-like haloes. Only quantities associated 
with the central galaxy and its main progenitors are plotted. Red lines show results from the Durham model, green lines show results 
from the Munich model and blue lines are for MORGANA. Solid and dashed lines in the second panels from the top show results for the 
tmrg = oo and tmrg = run, respectively. The haloes shown in this figure provide two representative examples for a halo with a rather 
quiet mass accretion history (left panels) and a halo with a larger number of massive mergers (right panels). 
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Figure 2. From loft to right: cold gas, hot gas, and gas cooling between two subsequent snapshots. Each line shows the average computed 
over 100 MW-like haloes. Solid and dashed lines in the left panels show results for the tmrg = oo and tmrg = run, respectively. Colour- 
coding is as in Fig. [l] 
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Figure 3. Cold gas for the tmrg = case (left panels), cold gas for the tmrg = oo case (middle panels) and hot gas (right panels) 
associated with the central galaxy for the MW sample. Top panels show results obtained by the Durham code against those from the 
Munich model. Bottom panels show results by the MORGANA code against those from the Munich model. The dashed line is the 
one-to-one relation and is plotted to guide the eye. 



Durham model, the cold gas fraction starts rising at z ~ 15, 
reaches values of ~ 0.09 at z ~ 10, and stays almost con- 
stant down to z — 0. The evolution of the hot gas fraction 
reflects that of the cold component at high redshift, with 
the Durham model predicting the largest hot gas fractions, 
due to a less efficient cooling with respect to the other two 
models. At z < 4, the Durham model is on average still 
slightly above the predictions from the Munich model and 
from MORGANA. On average, the Durham model predicts 
the highest cooling rates in the redshift interval 2 — 4. Predic- 



tions from the Munich model and from MORGANA are very 
similar, with a broad peak over the same redshift interval, 
and with a rapid decline at z < 2. 

Fig.[3]shows the amount of cold (left panels for tmrg = 
and middle panels for tmrg = oo) and hot gas (right pan- 
els) for all MW-like haloes in our sample, at z — 0. Top 
and bottom panels compare the Durham model and MOR- 
GANA with the Munich model, respectively. The agreement 
between the Munich model and MORGANA is very good for 
the predicted cold gas, particularly when zero merging times 
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are adopted, with only a very slight tendency for MORGANA 
to predict a smaller gas content with respect to the Munich 
model. When suppressing galaxy mergers {tmrg = oo), the 
agreement is still good, but the scatter is much larger and 
there is a systematic trend for larger cold gas amounts in 
MORGANA. The examples discussed above suggest that this 
might be due to a different treatment of the rapid accretion 
regime in these two models. The same systematic trend is 
visible for the hot gas content (right panel), which is always 
larger in MORGANA than in the Munich model. The trends 
are similar with the systematics being slightly stronger in 
the top panels, that compare prediction from the Durham 
model to those from the Munich model. This figure demon- 
strates that all three models used in this study provide very 
similar predictions at present, but that there are some resid- 
ual systematic trends at late times, and the evolution at high 
redshift is quite different (see Fig. [T] and Fig. [2}. 



5 SCUBA HALOES 

In order to investigate how the level of agreement between 
different models discussed in the previous section depends 
on halo mass, we have complemented our MW-like sample 
with a sample of SCUBA-like haloes. As discussed in Section 
(2] these have been selected as haloes with a number density 
of 10"'' at z ~ 2 and with a relatively massive descendant 
at 2; = (M200 = T.8 x lO" - 1.3 x 10^^ Mq). Our decision 
to use a SCUBA-like sample has been partially motivated by 
previous claims that MORGANA provides a good agreement 
with the observed redshift distribution and number counts 
of 850-/im sources because the adopted cooling model pre- 
dicts significantly larger coolin g rates w ith respec t to the 
'class ical' cooling model (Font anot et al. I [2007; Vi ola et al.l 
|2008| ). At this mass scale, we therefore expect significant 
differences between predictions from MORGANA and those 
from the Munich and Durham models, that both adopt dif- 
ferent implementations of the classical model. 

Fig. m shows the evolution of the cold and hot compo- 
nents (second and third panels from the top) normalized to 
the parent halo dark matter mass, and of the cooling rates 
for two examples from our SCUBA sample. As for Fig. [TJ 
these two examples have been chosen as representative for a 
halo with a rather quiet mass accretion history (left panels) , 
and one with a larger number of important mergers (right 
panels). The dark matter mass accretion histories of the 
haloes under consideration are shown in the top panels. Con- 
trary to what expected, predictions from MORGANA appear 
to be very close to those from the Munich model, while the 
Durham model deviates most from the other two, predicting 
systematically lower cold gas fractions at late times. The sys- 
tematics are stronger when satellite galaxies are allowed to 
survive as independent entities down to 2 = {tmrg = 00). 

The bottom panels of Fig. [3] show that the Durham 
model predicts much lower cooling rates than the Munich 
model and MORGANA, below z ~ 4. Over this redshift 
range, both the Munich model and MORGANA provide cool- 
ing rates as high as ~ 80 - 100 x 10^° MoGyr"^ in the two 
examples shown; values about twice (or more) as high are 
obtained in a number of the other haloes from the SCUBA 
sample. The behaviour predicted by the Munich model for 
these haloes appears to be somewhat noisier than that pre- 



dicted by MORGANA, with cooling rates that are not, how- 
ever, significantly lower. At the highest redshifts probed by 
the merger trees at disposal, the hot gas fraction predicted 
by the Munich model is lower than the corresponding value 
predicted by the Durham model and MORGANA. As noted 
in the previous section, this results from a very efficient cool- 
ing in the rapid accretion regime, as treated by the Munich 
model. The evolution of the hot gas fraction is then very 
similar in all three models, down to 2 ~ 2. At redshift lower 
than this, the Munich model and MORGANA still provide 
reasonably close results, while the Durham model is system- 
atically higher. 



The differences discussed above are clearly visible in 
Fig. O that shows the average evolution of the cold (left 
panel) and hot (middle panel) gas fractions associated with 
the central galaxy, and of the cooling rates, computed using 
all 100 haloes in the SCUBA sample. As in previous figures, 
solid and dashed lines in the left panel correspond to the 
tmrg = cxD aud tmrg ~ ruu, respectively. When satellite 
galaxies are merged immediately after their accretion, the 
cold gas fractions predicted by the Munich model and MOR- 
GANA are very similar at all redshifts, while the Durham 
model tends to predict larger cold gas fractions at high red- 
shift and lower fractions at z < 2. When satellite mergers 
are suppressed (i.e. satellite galaxies survive as independent 
entities, keeping the cold gas associated with them before 
accretion), the Durham model predicts significantly lower 
cold gas fractions than the Munich model and MORGANA. 
This is due to the significantly lower cooling rates predicted 
by the Durham model, as can be seen in the right panel of 
Fig. [5] This panel shows that, on average, the Munich model 
predicts the highest cooling rates on this mass scale, while 
predictions from MORGANA are intermediate between the 
Munich and Durham models. The middle panel of Fig. [5] 
shows the evolution of the hot gas fraction, normalized by 
the corresponding dark matter mass. It shows that, as noted 
for the MW-like haloes, the cooling efficiency at high red- 
shift is highest in the Munich model. This produces the rapid 
increase of cold gas visible in the left panel, and is due to a 
different treatment of the rapid cooling regime in this model. 
At z < 2, the Durham model predicts the largest hot gas 
fraction and the Munich model the lowest, with MORGANA 
again providing intermediate results. 



The amount of cold and hot gas for all haloes in our 
SCUBA sample at present are shown in Fig.[51 that compares 
predictions from Durham model and from MORGANA with 
results from the Munich model. From this figure, it is clear 
that there is a systematic trend for both the Durham model 
(top panel) and MORGANA (bottom panels) to predict lower 
cold gas fractions and higher hot gas fractions with respect 
to the Munich model. The disagreement is stronger for the 
Durham model than for MORGANA, and when galaxy merg- 
ers are suppressed. As shown above, this is due to systematic 
differences in cooling rates, which appear to be more signif- 
icant for this mass scale than for MW-like haloes (compare 
right panels in Fig. [5] and Fig. [2]). 
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Figure 4. Same as in Fig. [T] but for two SCUBA-like haloes. 
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6 NUMERICAL RESOLUTION AND 
SUBHALO SCHEMES 

In order to study how the results discussed in previous sec- 
tions are affected by numerical resolution, we have taken 
advantage of the mini-MSII. As explained in Section [2l this 
simulation has been run with the same initial conditions of 



the MS-II, but lower force and mass resolution (the same 
as for the MS). We have identified the same haloes used in 
our MW-sample and run each model on the corresponding 
merger trees. Haloes were matched across the two simula- 
tions by finding, at 2 = 0, all FOF groups in the mini-MS-II 
within a distance of lMpc/i~^ of the coordinates of the tar- 
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Figure 5. Same as in Fig. [J] but for the SCUBA sample. 
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Figure 6. Same as in Fig. [3] but for the SCUBA sample. The dashed line is the one-to-one relation and is plotted to guide the eye. 



get haloes from the MS-lfl For each halo in the original 
MW sample, the matched halo is found from this subset 
of haloes, as the one that minimizes the absolute value of 
Mvi, (mini-MSII) /Mvir(MS-II) - 1. We find that the matched 
halo lies within lOOkpcfe"^ of the target halo in more than 
90 per cent of the cases. In 76 per cent of the cases, the 
present virial masses differ by less than 10 per cent, while in 
95 per cent of the cases matched haloes have virial masses 
that differ by less than 20 per cent at present. 

Fig. [7] shows the mass accretion histories of two MW- 
like haloes as obtained from the MS-II (solid lines) and from 

* This search radius is a factor of ^ 7 smaller than the typical 
separation between lO^'^ ^<3- 



the mini-MSII (dashed lines). As for previous figures, mass 
accretion histories have been constructed by connecting each 
halo to its main progenitor (i.e. the most massive) at each 
node of the tree. At mini-MSII resolution, it is not possible 
to follow the evolution of these haloes past 2-^7 — 9, while 
the first progenitors of the MW-like haloes under consider- 
ation are identified at z ~ 14 — 15 at the resolution of the 
MS-II. Over the redshift range where haloes are identified 
in both simulations, the mass accretion histories extracted 
from the MS and mini-MSII are in quite good agreement. 
Most of the differences discussed below, should then be as- 
cribed to the ability to follow the evolution of the parent 
dark matter haloes up to higher redshift in the higher reso- 
lution simulation. 
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Figure 8. Evolution of the cold gas content (top panels), hot gas (middle panels) and the cooling rates (bottom lines) for one of the 
MW-like haloes used in this sample, at two different resolution levels. In each panel, solid lines show the evolution computed using the 
merger trees extracted from the MS-II, while dashed lines correspond to merger trees from mini-MSII. Left, middle and right panels 
show results from the Durham, Munich and MORGANA model, respectively. The mass accretion history of this halo is shown in the top 
panel of Fig. [7] 



Figs. [8] and [9] show the evolution of the cold and hot 
gas fractions (top and middle panels respectively) and the 
cooling rates (bottom panels) for the two examples whose 
mass accretion histories are shown in Fig. [T] In these figures, 
solid lines show the predictions from the Durham, Munich 
and MORGANA models obtained using the higher resolu- 
tion trees from the MS-II, while dashed lines show predic- 
tions from the lowest resolution simulations. For clarity, we 
have only shown the cold gas fractions obtained when merger 
times are set to infinity. Although the overall behaviour of 
cooling rates is similar in the two resolution runs, particu- 
larly at low redshift, it is clear from these figures that none 
of the models used in this study achieves a very good con- 
vergence. The bottom panels of Figs. [S] and [5] show that at 



halo appearance in the lowest resolution simulation, all three 
models generally predict a much larger cooling rate than 
obtained in the highest resolution simulation. The cooling 
rates adjust rapidly at approximately the levels predicted in 
the highest resolution simulation, but are always somewhat 
larger than the higher resolution predictions. As a conse- 
quence, cold gas fractions predicted using the lower resolu- 
tion trees are generally higher than obtained when using the 
higher resolution trees. 

In its standard implementation, the Munich model em- 
ployed in this study, runs on subhalo-based merger trees, 
rather than on FOF-based trees like those that we have used 
in previous sections. It is then interesting to ask how much 
results discussed above are affected by the use of a differ- 
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Figure 9. As in Fig. |8] The mass accretion history of this halo is shown in the bottom panel of Fig. [7] 



ent scheme for the construction of the input merger trees. 
To address this question, we have run the Munich model 
on all the subhalo trees corresponding to the FOF trees in 
both the MW and SCUBA samples discussed in Section [S] 
We remind the reader that the subhalo-base d merger trees 
for th e MS and MS-II were cons t ructed by ISpringel et al.l 
l|2005l) and iBovlan-Kolchin et al.l l|2009l ) as summarized in 
Section |21 and that they are publicly available. 



Figs. llOland llll show the evolution of the cold gas frac- 
tion (top panels), hot gas fraction (middle panels) and cool- 
ing rates (bottom panels) of the same MW and SCUBA-like 
examples shown in Figs. [1] and |4l In these figures, green lines 
correspond to results from the FOF-based trees (and are 
therefore equivalent to the results from the Munich model 
shown in Figs. [1] and |4| , while black lines show the corre- 
sponding results based on subhalo-based merger trees (red 



lines in Fig. [TT] will be discussed in Section [8] and can be 
ignored for now). 

These figures show that predictions obtained using the 
FOF-based trees are in quite nice agreement with those ob- 
tained using the subhalo-based trees. For MW-like haloes, 
there is some tendency for the FOF-based trees to provide 
lower cold and hot gas fractions with respect to results from 
subhalo-base trees. For the SCUBA haloes, the cold gas frac- 
tions predicted from the two schemes are very similar over 
all the redshift range where it is possible to trace the main 
progenitor of the haloes under consideration. At z < 2, the 
hot gas fraction predicted using the FOF-based trees is only 
slightly lower than results obtained using the standard (for 
the Munich model) subhalo scheme. 

The bottom panels of Figs. [TO] and [TT] show that the 
cooling rates obtained using the two difi^erent schemes are 
very close, particularly at high redshift. Some small differ- 
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Figure 10. From top to bottom: cold gas fraction, hot gas fraction, and gas cooling rate for the same MW-like haloes shown in Fig. [T] 



Solid and dashed lines in the top panels show results for the tmrg 
to subhalo-based and FOF-based merger trees, respectively. 



: oo and t„ 



: run, respectively. Black and green lines correspond 



ences are, however, visible. These arise from diflFerent halo 
merger times due to the subhalo scheme explicitly following 
haloes once they are accreted onto more massive systems. 
We recall that in the standard Munich model, a residual 
merging time (given by Eq. ^ is assigned to satellite galax- 
ies only when the dark matter substructures fall below the 
resolution limit of the simulation. We have kept this assump- 
tion in the examples shown here so that the tmrg = case 
corresponds to galaxies merging at the time their parent sub- 
haloes are tidally stripped at the resolution limit of the simu- 
lation, rather than instantaneously merging at the time their 
parent haloes become proper substructures. That is why the 
two dashed lines shown in the top panels of Figs. llOl and llll 
follow a different evolution. As we have explained in Section 
[S] in this model the amount of gas available for cooling is 
computed at the beginning of each snapshot by assuming 
baryon conservation (see Eq. [T]). Since then some cold gas is 
retained in distinct satellites for some time in the subhalo 
scheme, this leads to a different amount of gas available for 
cooling, and to the differences discussed above. It is interest- 



ing, however, that despite the delay due to the identification 
of dark matter substructures, the hot and cold gas fractions 
predicted by the two schemes, as well as the predicted cool- 
ing rates, do not differ dramatically. The differences between 
the two schemes become more important at lower redshift, 
but are in all cases smaller than those obtained from alterna- 
tive modelling schemes (see Sections |4] and [5]). This suggests 
that tidal stripping and truncation are very efficient in re- 
ducing the high redshift substructures below the resolution 
limit of the simulation, while dark matter subhaloes survive 
longer as independent entities at lower redshift. This is ex- 
pected, due to the increase of dynamical t imes at later times 
(see also Fig. 4 in lWeinmann et al.ll2009l 'l. 



7 MERGERS 

In this section, we compare the different implementations of 
galaxy mergers adopted in the three models used in this 
study. To this aim, we have identified all satellites that 
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Figure 11. From top to bottom: cold gas fraction, hot gas fraction, and gas cooling rate for the same SCUBA haloes shown in Fig. |4l 



Solid and dashed lines in the top panels show results for t„ 



CO and tr, 



0, respectively. Black and green lines correspond to 



subhalo based and FOF-based merger trees, while red lines correspond to the outputs from the Durham model, with darker lines showing 
results obtained assuming an isothermal profile for the hot gas. 



are present in each pair of modelijfl and that are assigned 
merger times that are lower than the Hubble time. As a 
reminder, merger times in the Durham model and in MOR- 
GANA are re-assigned after a new formation event, or halo 
major merger. For the comparison discussed below, we have 
always considered the last assignments in these two models. 

Fig. [12] compares the distributions of merging times 
in all three models used in this study. Top panels are for 
satellites of the MW-like haloes, while bottom panels are 
for SCUBA-like haloes. There is no significant difference be- 
tween the two samples, but a much larger number of satel- 
lites for the SCUBA haloes. Fig. [12] shows that a fraction 
of model satellites distribute along the one-to-one relation 



^ Due to a different treatment of the rapid cooling regime, there 
are haloes that host a galaxy in one model and not in another. 
These 'unpaired galaxies' are excluded in the analysis presented 
here. 



when comparing the merger times assigned in the Durham 
model with those obtained from the Munich model (left pan- 
els). The scatter is, however, very large, with a large number 
of satellites merging within relatively short times in the Mu- 
nich model and getting a much longer merger time in the 
Durham model, and viceversa. The correlation between the 
Munich model and MORGANA is much tighter (middle pan- 
els in Fig. I12|l . but there is a clear offset indicating that 
model satellites in MORGANA have merger times that are 
systematically lower than the corresponding times from the 
Munich model. Finally, when comparing the Durham model 
and MORGANA (right panels), all satellites fall below the 
one-to-one relation with a quite large scatter, and there is a 
large concentration of galaxies that merge within ~ 5 Gyr 
in the Durham model and within ~ 0.6 Gyr in MORGANA. 

At least part of the scatter in Fig. [12] is due to the 
random extraction of orbital parameters in the Durham 
and MORGANA models. To show a 'cleaner' comparison, we 
have re-run both models neglecting the orbital dependency. 
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Figure 12. Distributions of merging times in the three models used in this study. Top panels are for MW haloes while bottom panels 
are for the SCUBA haloes. Maps have been computed using only satellite galaxies present in both models and with merger time smaller 
than the Hubble time. 
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Figure 13. As in Fig. 1121 but neglecting the orbital dependencies in the merger time calculation for the Durham and MORGANA 
models (see text for details). 
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Figure 7. Mass accretion liistory of two of tlie haloes from the 
MW-like sample used in this study, at two different resolution 
levels. Solid lines show the histories constructed using the MS-II 
simulation, while dashed lines show the corresponding histories 
from the mini-MSII. 



In the Durham model, this has been achieved by setting 
©orbit = 1, while in MORGANA only circular orbits have 
been considered, with orbital energy set to 0.5 (see Section 
13. 3[) . Results are shown in Fig. 1131 As expected, the scatter 
is reduced in all panels. A larger fraction of model satellites 
distribute along the one-to-one relation in the left panels. 
There is, however, still a significant fraction of satellites that 
get longer merger times in the Durham model than in the 
Munich model. As discussed in Section 13.41 the dynamical 
friction formulations adopted in the Munich and Durham 
models differ by the assumptions made for the Coulomb 
logarithm, and for a different numerical pre-factor. All dif- 
ferences visible in the left panels of Fig. [T3] are due to these 
different assumptions, and to the re-assignment of merger 
times after each formation event (see below). 

There is still a strong correlation between the merger 
times assigned in MORGANA and in the Munich model, with 
an offset towards lower merger times in MORGANA. Neglect- 
ing the orbital dependency, the offset reduces at long merger 
times. The correlation between MORGANA and the Durham 
model is still quite bad, with a large fraction of satellites that 
merge within 0.6 Gyr in MORGANA and that are assigned 




Figure 14. Merger times (in units of dynamical times) for differ- 
ent mass ratios. Dashed, long-dashed and dot-dashed lines cor- 
respond to the standard assumptions adopted by the Munich, 
Durham and MORGANA models respectively. The dotted line 
has been obtained from the long-dashed line by changing the as- 
sumption for the Coulomb logarithm to A^f = 14- Mmain/j^^fsat- 
The predictions from MORGANA do not depend significantly on 
halo concentration. The thick solid line correspo nds to the fitting 
formula provided bv lBovlan-Kolchin et all ||2008| ). with no orbital 
dependency and for circular orbits. 



merger times that are about ten times longer in the Durham 
model. 

The results discussed above can be understood by com- 
paring Eqs.[3l[5l and[6l This is done in Fig. 1141 which shows 
the merger times obtained using these different formula- 
tions, for different mass ratios. In this figure, we have ne- 
glected the orbital dependencies and have assumed that halo 
concentration (this enters the equations adopted in MOR- 
GANA) varies as a function of mass according to Eq. 4 in 
iNeto et al.l(|2007l ). The figure shows clearly that the Munich 
and Durham implementations differ by a scaling factor but 
for Afaat/Mniain J? 0.1, duc to a different assumption about 
the Coulomb logarithm: the dotted line in Fig. 1 141 has been 
obtained from the long-dashed line by changing the assump- 
tion for the Coulomb logarithm to Adf = 1 -|- Mmain /A/sat, 
as in the Munich model. The re-calculation of merger times 
after each formation event or major merger apparently 
does not affect significantly the expected disagreement. It 
should be noted that in Figs. I12l and 1131 we have only used 
satellite galaxies with merger times lower than the Hub- 
ble times. These figures are therefore dominated by merg- 
ers with Msat/Mmain ^ 0.1. lu this regime, the MORGANA 
model is offset low with respect to the Munich model. It also 
predicts systematically lower merger times with respect to 
the Durham model but in this case the offset is not linear 
as a function of the mass ratio. In Fig. 1141 we have also 
shown the fitting formula provided bv lBovlan-Kolchin et al.l 
(200^, with no orbital dependency and for circular orbits. 
As noted earlier, the adoption of a simple fudge factor does 
not suffice to bring the formulation adopted in the Munich 



© 2010 RAS, MNRAS 000, [HgT] 



A SAM comparison - gas cooling and galaxy mergers 19 



model in agreement with the new fitting formula proposed 
by Boylan-Kolchin et al. 



8 SUMMARY AND DISCUSSION 

In this work we have compared results from three inde- 
pendently developed semi-analytic models, and we have fo- 
cused on alternative implementations of gas cooling and 
galaxy mergers. Our model comparison is carried out us- 
ing two large samples of identical merger trees, which allows 
us to compare results on a case-by-case basis and to fo- 
cus on differences due to different model assumptions and 
parametrizations. In particular, we have constructed two 
FOF-based merger tree samples. One is a set of 100 haloes 
from the MS-II with redshift zero dark matter halo masses 
similar to that of the Milky Way (MW-like haloes), while 
the other consists of 100 haloes from the MS with a num- 
ber density similar to that measured for SCUBA galaxies at 
z ^ 2 (SCUBA-like haloes). By using stripped-down versions 
of the models, we are able to avoid possible degeneracies and 
complications due to a different treatment of various physi- 
cal processes and to concentrate on the influence of specific 
assumptions and/or parametrizations. As explained above, 
we have chosen to include only the processes of gas cooling 
and galaxy mergers. In the following, we summarize briefly 
the results obtained and discuss them. 



8.1 Gas cooling 

Our results show that the different assumptions adopted 
about gas cooling lead to very similar results on mass scales 
similar to those of our own Galaxy, and to a generally good 
statistical agreement. Important differences arise, however, 
on larger mass scales. 

Two of the models used in this study (the Munich and 
Durham models) assu me variations of the cooling model 
originally proposed by IWhite &: FrenkI (|l991 ). The other 
model used here (MORGANA) adopts a more sophisti- 
cated model and, on the basis of previous published results 
j Viola et al.| [2008). was expected to provide systematically 
higher cooling rates for our SCUBA-like haloes. Contrary to 
this naive expectation, however, results from the Munich 
model and from MORGANA are very similar for this mass 
scale, while the Durham model gives systematically lower 
cooling rates. These results appear to be in contradiction 
with previously published tests, but the contradiction is only 
apparent as we discuss below. 

As explained in detail in Section [S] although the Mu- 
nich and Durham models are variations of the same cool- 
ing model, they differ in a number of details, in particu- 
lar for the assumption adopted for the hot gas distribu- 
tion. Red lines in Fig. [TT] show results from the standard 
Durham model (that assumes a /3 profile for the hot gas) 
and from a model that uses the same cooling implemen- 
tation but assumes an isothermal distribution for the hot 
gas (darker red lines), as in the Munich model. The figure 
shows clearly that by changing this assumption, the pre- 
dicted cooling rates are larger, bringing model results closer 
to (although they are still lower than) those obtained from 
the Munich model (green lines in Fig. Illf) . Some differences 



are still apparent, however: the Munich model predicts sys- 
tematically higher cooling rates than the Durham model, 
particularly at high redshift, where cooling is much more ef- 
ficient in this model than in both of the Durham implemen- 
tations considered here. These residual differences are due 
to a number of other different assumptions, in particular for 
the calculation of the cooling radius (see Section [3.41) . Sta- 
tistically, the differences between the models are relatively 
small, which reflects a general agreement in the underlying 
framework of these two models. 

Perhaps more surprising is the similarity between the 
results obtained from the Munich model and those from 
MORGANA, at b oth mass sca l es an alysed in this paper. As 
discussed earlier, IViola et al.l (|2008l 'l have claimed that the 
cooling model implemented in MORGANA predicts cooling 
rates that are significa ntly larger than t hose obtained using 
the classical model by I White fc FrenkI (|l99ll '). It should be 
noted, however, that their implementation of the classical 
model did not assume any special treatment for the 'rapid 
cooling regime' (as is instead done in the original work by 
White and Frenk and all subsequent variants of this model), 
and assumed that the hot gas distribution is described by a 
polytropic equation of state with index 7p = 1.15, which is 
similar to the /3 profile assumed in the Durham model. In ad- 
dition, results discussed in Viola et al. were obtained for iso- 
lated static haloes and it is not trivial to generalize them to 
the case of cosmological mass accretion histories, like those 
we have used here. Their conclusions is, however, valid when 
one adopts a gas distribution similar to that adopted in the 
Durham model used in this study (which is indeed similar 
to that they assumed in their implementation of the classi- 
cal model). Adopting a steeper gas profile, as in the Munich 
model used here, changes results significantly in some mass 
regimes, bringing them in very good agreement with those 
from the MORGANA model. 

It is important to realize that the differences highlighted 
above will have important consequences on predictions from 
galaxy formation models. For example, the Munich and 
MORGANA models will need to assume a stronger feedback 
than the Durham model, to counter-act excessive cooling 
at low redshift in relatively large haloes. Although all mod- 
els achieve this using very similar schemes (the AGN feed- 
back), the relative importance of these additional physical 
processes will be different in these three models. Predicting 
much lower cooling rates for massive haloes, the Durham 
model will have difficulties in providing large numbers of 
galaxies with elevated star formation at high redshift. In- 
deed, this model is able to reproduce the observed galaxy 
number counts at 850/im only by assuming a top-heavy ini- 
tial mass function (IMF) from the stars formed in bursts. 
This works because a top-heavy IMF has a much larger recy- 
cled gas fraction, which provides fuel for star formation. As 
illustrated above, MORGANA predicts much larger cooling 
rates than the Durham model and is indeed able to repro- 
duce the number counts of submilli miter sources (but th e 
brightest ones) with a standard IMF (jFontanot et al.llioOTi '). 
The Munich model used in this study predicts even larger 
cooling rates on average, but its predictions for the submil- 
limiter number counts have not been explored yet. 
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8.2 Galaxy mergers 

The three models used in this study adopt a different mod- 
elling for galaxy mergers. The Munich and Durham model 
assume variations of the classical dynamical friction formula. 
Results from these two models are somewhat correlated, 
but there is a large scatter and a large number of galax- 
ies get significantly longer merger times in one model than 
in the other. As explained in Section [T] this is mainly due 
to the adoption of a different numerical factor in front of 
the dynamical friction formula employed, and to a differ- 
ent assumption about the Coulomb logarithm. MORGANA 
uses formulae derived from numerical simulations and ana- 
lytic models that take into account dynamical friction, mass 
loss by tidal stripping, tidal di sruption of subhaloes, and 
tidal shocks (|Taffoni et al.ll2003l ). As shown above, over the 
range of mass-ratios that provide merger times lower than 
the Hubble time, the merger times obtained using these for- 
mulae are systematically lower (by a factor ~ 10 Tdyn) than 
those computed from the Munich model. Of the three mod- 
els used in this study, the Munich model uses merger times 
that are closer to the fitting formula recently proposed by 
iBovlan-Kolchin et al.l (2008), when neglecting any orbital 
dependence. 

The differences just discussed have important conse- 
quences for the stellar assembly history of massive galaxies, 
and for the formation and evolution of the brightest clus- 
ter galaxies and of the intra-group and intra-cluster light. 
MORGANA (and to some extent also the Durham model) 
will tend to assemble massive central galaxies earlier than 
the Munich model. To keep these galaxies red, these models 
will need to assume a somewhat stronger supernovae feed- 
back so as to make most of the mergers driving their late 
stellar assemble dry (i.e. avoid triggering late bursts that 
would rejuvenate the stellar population of these galaxies). 
A different balance between AGN cooling and tidal strip- 
ping of stars will also be required in these models to keep 
model predictions in agreement with observational results. 
These considerations are of course valid in the case all mod- 
els would use the same treatment of all other physical pro- 
cesses at play. We remind, however, that as mentioned in 
Section [T] these processes are usually treated in a different 
way complicating the comparison between different models. 



8.3 Numerical resolution and merger tree scheme 

Taking advantage of the mini-MSII, run using the same ini- 
tial conditions and volume as for the MS-11 but lower resolu- 
tions, we have analysed how the results discussed above vary 
as a function of numerical resolution. The level of agreement 
between the three models used in our study is not affected by 
numerical resolution. None of the models used in this paper, 
however, achieves a good numerical convergence, with all of 
them predicting moderately larger cooling rates in lower res- 
olution runs. 

On the other hand, results seem to be quite stable to al- 
ternative schemes for the construction of dark matter merger 
trees. In particular, we have compared results obtained using 
FOF-based trees with those obtained using subhalo-based 
trees, which represent the standard input of the Munich 
model used in this study. The small differences found by 
comparing results from these two schemes can be ascribed to 



the possibility of tracking the accreted haloes until they are 
stripped below the resolution limit of the simulation by tidal 
stripping and truncation. Interestingly, our results show that 
these processes act on relatively short time-scales, and that 
they are more efficient at higher redshift. 



8.4 What Next? 

One question that this study does not address is: what is the 
best way to model gas cooling and galaxy mergers? This is 
a very difficult question, particularly for the cooling model. 
It is clear that one needs to understand how the gas is dis- 
tributed in dark matter haloes, and how this distribution is 
affected by heating from supernovae and/or AGN feedback. 
Although hydrodynamical simulations of galaxy formation 
are becoming increasingly sophisticated, these physical pro- 
cesses still need to be included as 'sub-grid' physics, i.e. using 
prescriptions that are 'semi-analytical' in nature. As a con- 
sequence, published hydrodynamical simulations offer little 
indication of appropriate modelling of the hot gas distribu- 
tion and evolution. 

Our results have pointed out that modelling of the 
rapid cooling regime differ significantly in the implementa- 
tions discussed in this paper. Substantial numerical work 
has been focused recently on this mode of accret i on, a l- 
though it was was di scusse d as early as iBinnev (|l977h , 
iRees fc OstrikeJ j 19771 ). and IWhite fc FrenkI (|l99ll ). Inter- 
estingly, recent studies show that gas accretion during the 
'quasi-static regime' in hydrodynamical si mulations is sensi- 
tive to d ifferent implementat ions of SPH (|Keres et al.ll2009l . 
see also lYoshida et al] |2002| ) , while accretion rates in the 
rapid cooling regime are quite robust. One possible addi- 
tional concern in using numerical results to inform semi- 
analytic models is the poor performance of SPH codes in 
resolving and treating dynamical insta bilities developing a t 



resolvmg and treating dynamical msta bnities developing a t 
sharp interfaces in a multi-phase fiuid lAgertz et ahlboOTl ). 

The situation is somewhat better for the modelling 
of galaxy mergers. Since the merging process is predom- 
inantly driven by gravity, it can be studied using con- 
trolled numerical experime nts as done, for example, by 
iBovlan-Kolchin et al.l l|2008l ). As noted earlier, however, re- 
cent work has not yet converged on the dynamical friction 
formula appropriate for galaxy formation models. Further 
work in this area is therefore needed. 

Gas cooling and galaxy mergers are two basic ingredi- 
ents of any galaxy formation model that are relatively well 
understood. Also at this level, however, different assump- 
tions have to be made when implementing these processes. 
These give rise to non- negligible differences that can have 
important implications on the weight that needs to be given 
to additional physical processes (e.g. AGN feedback, tidal 
stripping of stars, etc.). This paper highlights specific areas 
where further work is needed in order to improve our galaxy 
formation models, with the ultimate goal of improving our 
understanding of the physical processes driving galaxy for- 
mation and evolution. 
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